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ABSTRACT 


A  simple,  microprocessor  oriented  algorithm  is  developed 
to  identify,  track  and  suppress  bending  mode  signals  from  a 
control  system's  rate  and  position  feedback  signals  using 
adaptive  digital  notch  filters.  The  algorithm  can  be  used 
to  suppress  bending  modes  having  center  frequencies  as  close 
as  one  octave  above  the  control  system  gain  cross-over  fre¬ 
quency  without  introduction  of  the  excessive  phase  loss 
associated  with  conventional  lowpass  filtering  techniques. 

A  third  order  model  of  the  trident  missile  autopilot  pitch 
attitude  control  loop  is  contaminated  with  two  dynamic, 
destabilizing  bending  modes  and  used  as  a  concept  demonstra¬ 
tion  model.  The  algorithm  is  demonstrated  by  stabilizing 
the  pitch  attitude  loop  in  the  presence  of  two  bending  modes 
with  unknown  gains,  damping,  center  frequencies  and  rates 
of  change  of  center  frequencies. 
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I .  INTRODUCTION 


The  purpose  of  this  thesis  is  to  develop  a  simple,  micro¬ 
processor-oriented  scheme  that  eliminates  body  bending  mode 
signals  in  a  missile  control  system  by  adaptive  digital  notch 
filtering  of  the  bending  mode  resonant  frequencies  from  the 
rate  and  attitude  gyro  feedback  signals.  Although  this  thesis 
uses  the  Trident  Fleet  Ballistic  Missile  as  a  conceptual 
demonstration  model,  the  process  is  applicable  to  many  control 
systems  having  unknown  and/or  time  varying  mechanical  reso¬ 
nances  . 

A.  THE  BODY  BENDING  PROBLEM 

Body  bending  is  a  composite  term  which  encompasses  the 
coupling  of  all  structure -borne  vibrations  from  various 
sources  into  the  missile  control  system  through  the  mechani¬ 
cal  sensors  such  as  rate  gyros  and  accelerometers.  The 
effect  of  the  body  bending  mode  signals  in  the  rate  and 
attitude  gyro  feedback  signals  is  destabilizing  since  if 
left  unattenuated,  the  missile  control  system  tries  to  com¬ 
pensate  for  the  sensed  body  bending  rate  through  the  thrust 
vector  control  (TVC)  system  thus  further  exciting  the 
bending  mode  and  the  inertial  or  rigid  body  rates  soon  be¬ 
come  insignificant  and  go  uncorrected. 

Some  of  the  sources  of  body  bending  are: 

1.  The  missile  acting  as  a  cantilevered  beam. 
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2.  Hydraulic  servo  resonance. 


5.  Structural  bending  at  the  attachment  points  for 
hydraulic  actuators. 

4.  Beam-like  resonance  of  sub-assemblies  such  as 
upper  stage  motors  and  payload. 

The  body  bending  problem  in  missiles  is  traditionally 
approached  in  one  of  three  ways: 

1.  The  mechanical  sensors  are  physicall)  located 
at  nulls  in  the  lower  frequency  bending  modes  and  the 
natural  bandwidth  of  the  control  system  and  lowpass  filters 
are  used  to  attenuate  the  higher  frequency  bending  modes. 
Reference  [1]  contains  a  detailed  description  of  the  analy¬ 
sis  and  testing  required  to  implement  this  approach  in  the 
control  system  of  the  XSM-65  missile.  The  popular  name,  if 
any,  of  the  XSM-65  missile  is  unknown. 

2.  The  mechanical  sensors  are  physically  located 
at  a  position  in  the  missile  where  the  deflections  due  to 

the  lower  frequency  body  bending  modes  are  in  phase  with  each 
other.  Phase  shifting  networks  are  used  to  ensure  that  the 
bending  mode  signals  are  fed  back  in  such  a  manner  that  the 
bending  modes  are  attenuated.  Lowpass  filters  are  still 
used  to  attenuate  the  higher  frequency  bending  modes. 
Reference  [2]  contains  a  general  description  of  how  this 
method  is  implemented  in  the  control  system  for  the  Saturn 
Missile. 
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Lowpass  or  bandstop  filters  are  used  in  the  feed- 


3 . 

back  control  loops  to  attenuate  the  body  bending  mode  fre¬ 
quencies  from  the  rate  and  attitude  gyro  feedback  signals. 

A  general  description  of  how  this  method  is  used  in  the 
Trident  Missile  Control  System  is  contained  in  Ref.  [3]. 

In  the  Trident  Missile  the  above  three  approaches  have 
many  undesirable  features.  Location  of  the  sensors  at  nulls 
or  where  the  body  bending  modes  are  in  phase  requires  a  very 
precise  knowledge  of  the  bending  mode  characteristics  and 
the  freedom  in  the  initial  design  to  locate  the  sensors  at 
the  anticipated  optimum  locations.  In  a  solid  fuel  missile 
such  as  the  Trident,  the  location  of  the  sensors  is  pre¬ 
determined  to  be  located  with  the  other  flight  control 
electronics  with  very  little  latitude  to  select  a  null  or 
position  where  the  bending  modes  are  in  phase.  The  second 
major  limitation  is  the  lack  of  precise  knowledge  of  the 
bending  mode  characteristics.  The  bending  mode  resonant 
frequencies,  gains,  damping  and  phase  are  a  function  of  many 
factors  including  time  into  the  flight,  attitude,  payload 
configuration,  skin  temperature,  and  thrust  level.  For 
example,  the  Trident  Missile  is  designed  to  compensate  for 
as  much  as  a  twenty-five  percent  variation  in  thrust  level 
due  to  slight  variations  in  the  solid  propellant  mixture  and 
casting  while  only  a  three  percent  thrust  variation  is 
acceptable  on  the  liquid  fueled  Saturn  Missile  [2].  Without 
a  precise  knowledge  of  the  parameters  controlling  the  bending 


mode  resonant  frequencies,  gains,  damping  and  phase,  the 
lowpass  and  bandstop  filters  must  be  designed  on  a  worst 
case  basis  and  as  a  result  reduce  the  control  system  band¬ 
width  and  phase  margins  much  more  than  if  a  narrow  notch 
filter  had  been  used  at  each  specific  bending  mode  frequency 
to  remove  the  bending  mode  component  from  the  rate  and 
position  gyro  output  signals. 

As  the  Trident  Missile  evolves  into  the  next  generation, 
the  anticipated  bending  mode  frequencies  will  become  much 
lower  and  approach-  -he  natural  frequencies  of  the  flight 
control  system.  The  use  of  conventional  lowpass  filters 
under  these  conditions  would  produce  so  much  phase  loss  that 
the  missile  would  become  unstable. 

B.  THE  PROPOSED  SOLUTION 

The  proposed  solution  to  the  bending  mode  problem  is  an 
extension  of  preliminary  work  done  by  Lockheed  Missiles  and 
Space  Company  [3]. 

If  the  bending  mode  resonant  frequencies  can  be  deter¬ 
mined  in  near  real  time  during  the  flight  and  tracked  as 
their  controlling  parameters  change,  then  adaptive  notch 
digital  filters  in  the  rate  and  position  gyro  signal  feed¬ 
back  paths  can  be  implemented  to  attenuate  a  selected  number 
of  bending  mode  signals  that  are  contaminating  the  desired 
rigid  body  rate  and  attitude  signals. 


The  difficulty  in  the  implementation  of  such  a  solution 
lies  in  the  in-flight  determination  of  the  bending  mode 
frequencies  rapidly  and  with  sufficient  accuracy.  The  rapid 
determination  of  the  bending  mode  frequencies  is  essential 
in  order  to  place  the  notch  filter  at  the  proper  bending 
mode  frequency  before  the  bending  mode  signals  can  degrade 
the  missile  attitude  and  rates  to  a  point  from  which  it 
cannot  recover  within  the  allowable  flight  trajectory 
envelope.  The  allowable  time  for  the  initial  determination 
of  the  bending  mode  frequencies  is  a  function  of  the  initial 
frequency  and  bandwidth  of  the  individual  notch  filters. 

In  general  the  time  available  to  eliminate  the  effects  of 
the  bending  modes  before  the  Trident  Missile  deviates  from 
an  acceptable  flight  trajectory  is  one-half  second.  The 
requirement  for  accuracy  is  essential  to  allow  the  minimum 
bandwidth  notch  filters  to  be  used.  As  successive  notch 
filters  are  added  to  the  control  system  to  attenuate  the 
bending  mode  signals,  an  increasing  amount  of  phase  loss 
will  be  experienced  at  the  frequencies  used  by  the  control 
system  resulting  in  a  decrease  in  the  already  minimal  phase 
margin  for  overall  missile  stability.  If  very  narrow  notch 
filters  are  used  to  minimize  the  phase  loss  at  the  control 
system  frequencies,  then  the  center  frequency  of  the  notch 
must  be  placed  much  more  accurately  to  be  effective. 
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Four  methods  were  investigated  as  to  their  potential  to 
be  used  for  inflight  bending  mode  frequency  identification 
and  tracking: 

1.  Reference  [3]  summarizes  the  use  of  Fourier 
integral  filters  to  measure  the  power  spectral  density  at 
integer  frequencies.  Second  order  curve  fitting  techniques 
were  then  used  to  estimate  the  bending  mode  frequencies 
from  the  output  of  the  filters.  The  major  drawback  of  this 
technique  is  that  it  required  at  least  one -half  second  of 
integration  time  for  each  repeated  frequency  determination. 
Reference  [3]  did  not  address  the  ability  of  this  technique 
to  accurately  track  a  dynamic  bending  mode,  although  the 
data  indicated  that  a  bending  mode  that  was  changing  as 
little  as  0.2  hertz  per  second  would  cause  sufficient  dis¬ 
tortion  of  the  frequency  determination  as  to  render  the 
technique  ineffective. 

2.  The  author  conducted  a  feasibility  study  of 
using  a  Fast  Fourier  Transform  (FFT)  algorithm  to  determine 
the  spectrum  of  the  rate  gyro  output  signal  in  order  to 
determine  the  bending  mode  frequencies.  The  method  was 
disregarded  for  two  reasons.  The  rate  gyro  signal  would 
have  to  be  sampled  for  ten  seconds  to  obtain  a  minimum 
frequency  resolution  of  0.1  hertz  and,  the  computational 
capability  required  to  execute  the  FFT  algorithm  was  felt 
to  be  excessive. 
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The  feasibility  of  using  a  phaselock  loop  to 


m. 


3 . 

lock-on  and  track  a  bending  mode  frequency  was  investigated 
by  the  author.  Simulation  of  a  phaselock  loop  with  an  eight- 
hertz  center  frequency  and  a  three-hertz  acquisition  range 
revealed  that  while  the  tracking  characteristics  of  the 
phaselock  loop  were  very  good,  the  ability  to  achieve  lock 
with  a  random  phase  signal  at  the  extremes  of  its  acquisition 
range  was  unreliable  if  the  signal  was  contaminated  with  fre¬ 
quency  components  from  the  plant  and  a  second  bending  mode 
at  anticiapted  signal  levels.  The  technique  was  not  pursued 
further  but  with  the  addition  of  prefilters  on  the  phase- 
lock  loop  input  this  technique  has  the  potential  to  be  more 
effective  than  the  technique  used  in  this  thesis.  Consi¬ 
derable  refinement  in  the  simulation  of  the  phaselock  loop 
nonlinearities  would  be  reqiiired. 

4.  The  technique  used  in  this  thesis  to  identify 
and  track  the  bending  mode  frequencies  uses  a  microprocessor 
to  compute  the  frequency  of  each  bending  mode  by  use  of 
digital  bandpass  filters  and  linear  estimation  of  the  zero 
crossing  times.  The  use  of  discrete  signal  processing  by 
the  microprocessor  will  not  require  additional  sampling 
circuitry  since  the  missile  attitude  and  rate  gyro  signals 
are  currently  sampled  and  put  in  digital  form  in  the  Trident 
Missile  for  use  in  the  digital  autopilot.  As  shown  in 
Figure  1-1,  the  microprocessor  uses  the  sampled  rate  gyro 
data  by  passing  it  through  a  set  of  parallel  digital 
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Figure  1-1  Microprocessor  Algorithm  Modules 
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bandpass  filters  in  order  to  decompose  the  signal  into  fre¬ 
quency  bands  that  are  only  wide  enough  to  contain  one 
bending  mode  frequency.  The  output  of  each  of  the  bandpass 
digital  filters  is  examined  after  each  sample  to  see  if  the 
amplitude  of  the  signal  in  that  band  has  changed  signs.  If 
the  amplitude  has  changed  signs  then  the  time  of  the  zero 
crossing  is  linearly  estimated  based  on  the  amplitude  of  the 
signal  before  and  after  the  zero  crossing  and  the  sampling 
period.  The  half  period  detected  is  converted  to  a  fre¬ 
quency  and  averaged  with  the  frequency  that  the  digital 
notch  filter  devoted  to  that  band  is  currently  tracking. 

The  average  of  the  current  and  computed  frequency  is  then 
used  to  compute  new  coefficients  such  that  the  digital  notch 
filter  will  be  centered  at  the  new  average  frequency.  The 
details  of  the  algorithm  and  design  are  discussed  in  Chapter 
Three  and  the  concept  is  demonstrated  by  simulation  in 
Chapter  Four. 


II.  MODEL  DEVELOPMENT 


In  order  to  verify  the  ability  of  the  tracking  notch 
filter  concept  to  stabilize  the  missile  in  the  presence  of 
bending  mode  feedback,  a  model  of  the  system  was  developed 
to  be  used  as  the  basis  for  simulation. 


A.  THE  MODEL  OF  THE  PLANT 

Figure  2-1  is  a  block  diagram  representation  of  the  re¬ 
duced  order  pitch  plane  attitude  control  loop  of  the  Trident 
I  Missile  autopilot.  The  model  was  provided  by  Lockheed 
Missile  and  Space  Company  by  reference  [3] ,  and  was  modified 
only  to  the  extent  that  the  signs  at  the  summers  for  and 
were  changed  to  allow  all  the  model  coefficients  to  be 
considered  positive  quantities.  Table  2-1  shows  the  time 
varying  values  for  the  parameters  used  in  this  model.  The 
model  shown  in  Figure  2-1  represents  the  rigid  body  dyna¬ 
mics  of  the  missile  and  is  referred  to  as  the  plant.  When 
bending  modes  are  added  to  the  rigid  body  model  the  resulting 
composite  model  is  referred  to  as  the  system. 

The  state  space  representation  of  the  plant  for  any  time 


is  : 


X1  =  X2 


x,  -  -K9MrX1-KrQKeM5X2-MaX5-K0M50c 


IS 


TABLE  2-1 


MISSILE  PARAMETER 

.  VALUES  AT 

VARIOUS  TIMES 

INTO  THE  FLIGHT 

PARAMETER 

T=4 0  SEC. 

T=  50  SEC. 

M 

a 

13.625 

28.625 

17.5 

22.51 

24.17 

32.0 

7 

“a 

88.87 

180.0 

140.0 

“5 

119.19 

113.75 

137 . 5 

V 

1750 

2500 

3400 

*r* 

0.06 

0.06 

0.06 

^GA 

0.06 

0.06 

0.06 

CD 

1.5 

1.48 

0.946 

K  a 
r9 

0.2 

0.2 

0.2 

NOTE:  The  units  of  each  parameter  are  listed  on  the  list 

of  symbols.  The  signs  of  the  parameters  are  all 
positive  to  be  consistent  with  the  structure  of  the 
model  used. 


t 

V 

\ 


x3  =*K9“5xi+C1'KreK9^s^x2'^a:c3'K0i:o< 


V 


V 


V 


where : 

is  the  missile  pitch  angle. 

x7  is  the  missile  pitch  angle  rate. 

x,  is  the  angle  of  attack. 
j  ° 


Using  matrix  notation  the  plant  is  then  described  by 

X  =  AX  +  B0 
~  ~~  ~  c 

G  =  CX 

where 


0 

1 

0 

-Vs 

-KreK9M5 

-M 

a 

■hh 

(1'Kr9K9^^ 

-  7 

V 

V 

V 

[°  -k9m5 

T 

K  Z  1 

9^5  J 

V 

[1  0 

0] 
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The  transfer  function  relating  the  missile  pitch  angle 
9  to  the  commanded  pitch  angle  9 is: 

6(s) 

-  =  H  ( s) 

9cCs) 


H(s)  =  -3 


V  V' 


s*'+(5i+KpKr9Ms)sZ+[Mw  +  KvftKa(M^-M,Z,)+KaM<r]s+KaCMs^>M 


a  r9  9  5-ct  “a-5 ;  ‘'9“5J  *  “9  ^‘5^  “a 
V  V  V 


Cl  J 


The  transfer  function  can  be  factored  with  a  root  finder 
and  used  to  plot  the  poles  and  zeros  of  the  transfer  function 
in  order  to  provide  an  intuitive  feeling  as  to  the  response 
of  the  system  not  otherwise  apparent  with  the  state-space 
representation.  The  state-space  representation  of  the  sys¬ 
tem  is  used  primarily  in  the  simulation  programs  to  obtain 
the  time-domain  solution  to  the  differential  equations  of 
motion. 

B.  THE  CHARACTERISTICS  OF  THE  PLANT 

The  model  of  the  plant  given  in  Figure  2-1  is  a  modelvof 
the  rigid  body  dynamics  including  two  Newtonian  integrators 
and  the  positive  feedback  of  the  angle  of  attack,  x^.  The 


„  v 
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output  pitch  angle  9  is  opposite  in  sign  to  the  commanded 


pitch  9  .  This  reversal  is  compensated  for  in  the  outer 
guidance  loop  not  included  in  this  model.  The  sign  of  the 
forward  gains  and  Z^.  are  normally  specified  as  negative 

quantities  which  results  in  a  net  negative  feedback.  For 
computational  convenience  the  gains  are  treated  as  positive 
quantities  and  the  signs  at  the  summers  have  been  changed. 

The  sign  of  the  rate  and  position  feedback  summers  remains 
positive  to  clearly  indicate  the  nature  of  the  positive 
feedback  of  the  bending  modes  from  the  pitch  rate  and  pitch 
angle  gyros. 

The  model  of  the  plant  was  designed  to  simulate  the 
transient  behavior  of  the  system  and  not  the  steady  state 
behavior.  As  a  result  the  simulation  studies  are  only  valid 
for  about  two  seconds  of  transient  behavior. 

To  establish  a  baseline  or  nominal  performance  the  plant 
was  simulated  for  a  ramp  and  a  step  input  with  system  para¬ 
meter  values  fixed  at  the  time  =  40  second  values.  The 
reasoning  behind  the  selection  of  inputs  and  parameter  values 
is  discussed  in  Chapter  Four. 

Selecting  the  appropriate  values  for  the  system  para¬ 
meters  from  Table  2-1  for  time  =  40  seconds,  the  transfer 
function,  given  by  equation  (1)  becomes: 

-55.77(3+0.01779) 

Hfs)  - - (2) 

($+0.0098) (S+3.603+/-J  7.17) 
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The  poles  and  zeros  of  the  transfer  function  are  plotted 


in  Figure  2-2. 

An  examination  of  the  poles  and  zeros  of  the  transfer 
function  indicates  the  system  is  stable  with  a  pair  of  com¬ 
plex  poles  having  a  natural  frequency  =  8.02  RAD/SEC  and 
a  coefficient  of  damping  £  =  0.449.  There  is  also  a  very 
slow  pole  at  s  =  -0.0098  and  a  zero  at  s  =  -0.01779. 

To  establish  the  baseline  or  nominal  performance  of  the 
plant  without  bending  modes,  the  plant  was  simulated  using 
the  IBM  DSL/360  language.  The  values  of  the  plant  para¬ 
meters  from  Table  2-1  for  time  =  40  seconds  are  used  through¬ 
out  this  thesis.  A  detailed  discussion  of  the  simulation 
programs  and  the  selection  of  the  plant  parameter  values  is 
contained  in  Chapter  Four.  Figure  2-3  illustrates  the  plant 
response  to  a  step  input  pitch  command  of  0.1  radians. 

Figure  2-4  illustrates  the  response  of  the  plant  to  a  ramp 
input  of  0.05  radians/sec. 

Note:  It  will  be  the  convention  in  this  thesis  when 

presenting  computer  generated  output  in  the  form  of  a  figure 
to  include  the  computer  program  that  was  used  to  generate 
that  figure  in  the  computer  program  appendix  properly 
labelled  to  indicate  the  figure  number.  If  reference  is 
made  to  the  program  it  will  be  referred  to  as  program 
(figure  number) . 
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Time,  Sec. 


Figure 


Sys  tem : 
Input : 


Third  order  missile  model  without 
bending  modes 

0.1  Radian  step  pitch  command 


2-3  Missile  response  for  a  step  command 


Although  all  simulations  in  this  thesis  are  performed 
using  the  third-order  model  of  the  plant,  it  is  possible  to 
reduce  the  model  to  a  second-order  model  without  introducing 
significant  error.  The  second-order  model  provides  a  con¬ 
venient  preliminary  design  model  since  the  characteristics 
of  second-order  systems  are  wrell  established.  The  second- 
order  model  is  derived  by  analysis  of  the  time  domain  res¬ 
ponse  of  the  third-order  system  and  eliminating  the  effects 
of  the  slow  pole  and  zero.  The  validity  of  the  second- 
order  model  is  established  by  comparison  with  the  third-order 
model  simulation  results. 

To  obtain  the  time  domain  response  to  a  0.1  radian  step 
input  the  transfer  function  of  equation  (2)  will  be  multi¬ 
plied  by  0.1/s  and  the  resulting  function  expanded  in  a 
partial  fraction  expansion. 


0.1  -5. 577 (s+0. 01779) 

9  (s)  =  -  H  (s )  =  - 

s  s(s+0.0098) (s+3.608+/-J  7.17) 


which  yields: 


-  0 . 0098t  -  3 . 608t 

0 (t )  «  -0.1  +  0.0453  e  +  0.06218e  COS (7 . 1 7t - 0 . 453) 

at  t  =  2.0  sec. 

9  (2)  =  -0.1  +  0.0453  (0.98)  +  0 . 062 1 8  (0 . 000 72 ) COS (1 4 . 34  -  0 . 4 55) 


=  -0.552 
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Examining  the  time  domain  response  one  finds  that  the 
influence  of  the  slow  pole  is  so  slow  that  one  may  consider 
the  response  of  the  system  due  to  that  pole  to  be  essentially 
a  constant  for  the  two  second  duration  of  the  simulation  with 
less  than  two  percent  error.  Treating  the  term  due  to  the 
slow  pose  as  a  constant  reduces  the  time-domain  solution  to: 

-3. 608t 

9(t)  =  -0.0547  +  0 . 06 21 8e  COS(7.17t  -  0.453)  (3) 

The  system  now  can  be  considered  as  second-order  and  the 
response  can  be  very  closely  approximated  from  tabulations 
of  the  characteristics  of  second  order  systems  [4] . 

The  transform  of  equation  (3)  is 

L[9(t)]-e(s)-  ^-°-5-4-  -  i0.0.5..589),(s.l7.098) 

S  s  +  7.226s  +  64.4265 


Placing  all  terms  over  a  common  denominator: 


9(S)  = 


-0.0547  s2-0.3952  s-3.524  +  0,05589  s2  +0.5967  s 
s  (s"+7.226  s+64.4265) 


Cancelling  the  essentially  equal  s"  and  s  terms  in  the 
numerator  and  factoring  out  the  transform  of  the  0.1  radian 
step  input  (0.1/s)  yields  the  following  reduced  order  trans¬ 
fer  function  relating  the  missile  pitch  angle  9  to  the 
commanded  pitch  angle  3  . 


’  HrCsl  * 


-35.24 


s  +  7.226  s  +  64.4265 


_  s 


The  steady  state  gain  of  the  system  is: 


s  H  (s ) 
r  ^  1 

LIM  -  =-0.547 

S+0  s 


The  time  of  peak  overshoot  will  be 

t  =  - — -  =  0.438  sec 

P 

W 

n 

Using  the  approximation  for  the  settling  time 

.  4.0 

t  =  -  =  1.11  sec 

S  cw 

n 

The  validity  of  the  second-order  model  for  a  ramp  input 

is  established  by  comparison  of  the  value  of  the  missile 

pitch  angle  determined  analytically  from  the  second-order 

transfer  function  at  time  =  2  seconds  with  the  value  resulting 

from  the  simulation  of  the  third-order  plant. 

Given  a  transfer  function  in  the  form: 

B„  +  B.s  +  B7s2  +  .  +  Bsm 

0  12  m 

G (s)  =  - 

A0  *  A lS  *  A2s2  *  .  *  Ansn 
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The  steady  state  time  domain  response  of  the  output  to  a 
ramp  input,  u(t)  =  at,  can  be  determined  from  the  following 
formula  [ 5] : 

B0  (A0B1  '  Vi> 

y_ ct)  *  —  at  ♦  - 

s  A  2 

0  Ao 

Using  the  second-order  transfer  function  determined  in 

equation  (3a)  one  finds  that  the  output  at  time=2  for  a  ramp 

input,  u(t)  =  O.OSt  is: 

-  35 . 24 (0 . 05) t  0.05  (35.24) (7.226) 

9  (t)  =  -  +  - » -  =  -0.05242 

s  64.4265  (64 .4265; 

Table  2-2  summarizes  the  results  obtained  from  the  second- 
order  model  with  those  resulting  from  simulation  of  the  third- 
order  model.  Again,  it  must  be  noted  that  the  simulation 
studies  contained  in  this  thesis  were  performed  using  the 
third-order  model  of  the  plant.  The  second-order  model  was 
derived  only  as  a  convenient  preliminary  design  tool  for  any 
subsequent  studies. 
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TABLE  2-2 

COMPARISON  OF  THIRD-ORDER  AND  SECOND-ORDER  MODEL  RESPONSES 


CHARACTERISTIC 

SECOND- 

ORDER 

MODEL 

THIRD- 

ORDER 

SIMULATION 

MAX I MUM 

OVERSHOOT 

1.21 

1.21 

TIME  AT  MAX. 

OVERSHOOT 

0.44  SEC 

0.44  SEC 

STEADY  STATE 

GAIN 

-0.5625 

-0.547 

SETTLING 

TIME  (2%) 

1.11  SEC 

1.11  SEC 

RAMP  OUTPUT 

AT  T=  2  SEC 

-0.0524 

-0.0526 
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C.  THE  MODEL  OF  THE  BENDING  MODES 


Reference  [3]  provided  the  model  of  the  plant  with  one 
bending  mode  shown  in  Figure  2-5.  The  bending  mode  attitude 
signal  x^  and  bending  mode  rate  signal  are  added  to  the 
rigid  body  attitude  signal  x^  and  rigid  body  rate  signal 
respectively.  The  composite  attitude  and  rate  signals  are 
each  fed  back  to  the  system  input  in  a  positive  sense  and 
the  overall  system  becomes  highly  unstable. 

Reference  [2]  states  that  the  bending  mode  contamination 
of  the  rigid  body  rate  signal  has  a  much  greater  destabilizing 
effect  than  does  the  bending  mode  contamination  of  the  rigid 
body  attitude  signal.  The  bending  mode  model  is  therefore 
based  on  modeling  the  bending  mode  rate  and  the  bending  mode 
attitude  signal  is  formed  by  integrating  the  rate  signal. 

Reference  [6]  contains  the  bending  mode  frequency  res¬ 
ponse  for  a  Poseidon  Missile  when  the  bending  modes  were 
excited  by  slewing  the  first  stage  TVC  nozzles.  Although 
the  data  was  measured  on  the  Poseidon  Missile  with  various 
flight  configurations,  the  results  are  consistent  with  the 
general  distribution  of  bending  modes  shown  in  Figure  3.3-1 
of  Ref.  [2],  The  distribution  of  the  bending  modes  indicates 
that  the  bandwidth  of  the  bending  modes  increases  as  the 
center  frequency  of  the  bending  mode  increases.  The  distri¬ 
bution  of  the  bending  modes  used  in  this  thesis  is  shown  in 
Figure  2-6.  The  low  frequency  bending  mode  may  be  centered 
anywhere  in  a  band  from  5  to  8  hertz  and  the  high  frequency 
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Figure  2-5  Model  of  the  Plant  with 
One  Bending  Mode 
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CONTROL  FREQUENCY  RANGE 


Figure  2<-6  Dynamic  Bending  Mode  Distribution 


bending  mode  may  be  centered  anywhere  in  the  band  from  8  to 
15  hertz.  The  high  frequency  bending  mode  has  approximately 
twice  the  bandwidth  of  the  low  frequency  bending  mode.  The 
simulation  studies  contained  in  Chapter  Four  demonstrate  the 
ability  of  the  adaptive  notch  filtering  scheme  to  identify, 
track  and  filter  bending  modes  that  have  time  varying  gains, 
bandwidths  and  center  frequencies. 

The  bending  mode  rates  were  modeled  as  the  outputs  of 
bandpass  filters  with  a  very  high  Q.  The  transfer  function 
relating  the  output  Y(s)  to  the  input  U(s)  for  a  general 
second-order  bandpass  filter  is: 

Bs 

GCs)  =  - 

s“+Bs+W“ 

c 

where : 

Wc  =  the  center  frequency  of  the  filter 

B  =  the  5  dB  bandwidth  of  the  filter 

The  block  diagram  of  G(s)  is  shown  in  Figure  2-7.  Y(s) 
represents  the  bending  mode  rate  output  and  since  the  bending 
mode  attitude  output  is  the  integral  of  the  bending  mode  rate, 
an  output  of  the  bending  mode  attitude  is  also  available  from 
the  rate  model  as  indicated  in  Figure  2-8  where  Y^fs)  is  the 
bending  mode  rate  output  and  Y^Cs)  is  the  bending  mode  atti¬ 
tude  output. 
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Inspection  of  Figure  2-8  reveals  that  modeling  the 
bending  mode  rate  as  a  bandpass  filter  is  equivalent  to  the 
bending  mode  model  shown  in  Figure  2-5. 

Figure  2-9  shows  the  system  model  consisting  of  the  plant 
and  two  bending  modes  that  are  used  for  the  simulation 
studies  in  this  thesis. 

D.  CHARACTERISTICS  OF  THE  PLANT  WITH  TWO  BENDING  MODES 

Preliminary  values  for  the  bending  mode  parameters  of  gain, 
center  frequency  and  bandwidth  were  selected  for  the  purposes 
of  analysis  and  to  develop  some  insight  as  to  the  elements 
controlling  the  behavior  of  the  system.  The  bending  mode 
parameters  are  varied  during  the  simulation  studies  presen¬ 
ted  in  Chapter  Four  but  the  detailed  analysis  is  not  repeated 
for  each  parameter  study. 

The  bending  mode  parameters  selected  for  analysis  of  the 
system  are: 

Gh  =  20.0 
Bh  =1.10  RAD/ SEC 

W^  =  62.83  RAD/SEC  (10  Hz  center  frequency) 

G^  =  240.0 
Bl  =0.1  RAD/SEC 

W  =  37.7  RAD/SEC  (6  Hz  center  frequency) 

The  frequency  response  of  the  transfer  function  relating 
the  missile  pitch  angle  9(s)  to  the  input  pitch  command  0c(s) 
with  the  rate  and  attitude  feedback  loops  open  was  determined 
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for  various  configurations  of  the  plant  and  bending  inodes 
shown  in  Figure  2-9.  The  plant  parameters  from  Table  2-1 
for  time  =  40  seconds  were  used. 

Figure  2-10  shows  the  open  loop  frequency  response  of  the 
plant  with  no  bending  modes. 

Figure  2-11  shows  the  open  loop  frequency  response  of  the 
6.0  Hz  bending  mode  alone. 

Figure  2-12  shows  the  open  loop  frequency  response  of  the 
plant  and  the  6.0  Hz  bending  mode. 

Figure  2-13  shows  the  open  loop  frequency  response  of  the 
10.0  Hz  bending  mode  alone. 

Figure  2-14  shows  the  open  loop  frequency  response  of  the 
plant  and  the  10.0  Hz  bending  mode. 

Figure  2-15  shows  the  open  loop  frequency  response  of  the 
plant  with  both  the  6.0  Hz  and  the  10.0  Hz  bending  modes. 

The  matrix  representation  of  the  state  equations  derived 
from  the  closed  loop  system  shown  in  Figure  2-9  is: 

« 

x  =  Ax  +  B6 
-  c 

9  =  Cx 

where : 

is  the  rigid  body  pitch  angle  of  the  missile 

x?  is  the  rigid  body  pitch  rate  of  the  missile 

x,  is  the  missile  angle  of  attack 

x.  is  the  pitch  angle  resulting  from  the  lower  frequency 
bending  mode 


_  v 
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Frequency,  He  (linear  scale) 


Figure  2-12  Open  Loop  Frequency  Response  of 
the  Plant  and  Low  Frequency 
Bending  Mode 
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Frequency,  Hz  (linear  scale) 


Open  Loop  Frequency  Response  of 
the  Plant  and  High  Frequency 
Bending  Mode 


Xp  is  the  pitch  rate  resulting  from  the  lower  frequency 
bending  mode 

x.  is  the  pitch  angle  resulting  from  the  higher  fre¬ 
quency  bending  mode 

x._  is  the  pitch  rate  resulting  from  the  higher  frequency 
bending  mode 

0c  is  the  commanded  pitch  angle 

S  is  the  missile  pitch  angle 


The  A  matrix  is : 


All 

A12 

A13 

A14 

A15 

A16 

A17 

A21 

A2  2 

A23 

A24 

A2  5 

A26 

A27 

A51 

A3  2 

A3  3 

A34 

A3  5 

A3  6 

A3  7 

A41 

A42 

A4  3 

A44 

A4  5 

A46 

A47 

A51 

A5  2 

A53 

AS4 

A55 

A56 

A57 

A61 

A62 

A63 

A64 

A65 

A66 

A67 

A71 

A72 

A73 

A74 

A75 

A76 

All 

The 

B  matrix  is: 

[B1 

B2  B3 

B4 

B5  B6 

B7]T 

The 

C  matrix  is : 

[Cl 

C2  C3 

C4 

C5  C6 

C  7  ] 
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where  : 


All  = 

A12  = 

A13  = 
A14  = 
A15  = 
A16  = 
A1  7  = 

A4 1  = 

A42  = 

A43  = 
A44  = 

A4  5  = 

A46  = 
A4  7  = 

A71  = 

A72  = 

A73  = 


0.0 

A21 

- 

-K9M3- 

A31 

1.0 

A22 

= 

-Kr9K9Md- 

A3  2 

0.0 

A23 

M 

A33 

a 

0.0 

A24 

a 

'^GA^d^S 

A34 

0.0 

A2S 

a 

’^rGZr9Z9M5 

A3  5 

0.0 

A26 

a 

‘^gaV’s 

A36 

o 

o 

A27 

"^rGZr9  Z6M5 

A37 

K  Zi 
0— 

1  -  K  „K„Z<$. 


re  0 


T 


a 

V 


i  *1 

7 

~  '  5 

■vrG-r9fae V 

.*  :** 
^GA“'9lT 


rG“r8  “eT" 


0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 


AS1  * 

A5  2  * 

A53  * 
A54  « 

A55  * 

AS  6  = 
AS7  = 


lbl 

A61  =  0.0 

GA 

K, 

R  If  6 

L  L  r9*GA 

A62  =  0.0 

.0 

A64  =  0.0 

LBLK9-WL2 

A65  =  0.0 

Gt  B  ip  rK  oT-— Bt  A66  =  0.0 
L  LyrG  L 

G.B.K.  A66  =  0.0 

*rC 

GTBrK  rtK„-~  A67  =  1.0 


L  L  r91ij) 


GA 


C’hBh*ga 

K 

B1  =  0.0 

Cl  =  1.0 

Wre*^ 

B2  =  K9M5 

7 

C2  =  0.0 

0.0 

=  -  Ke^ 

C3  =  0.0 
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A74  = 

GHBHK9 

<|> 

B4  = 

0.0 

G.  Bt  ,  9 

L  LuGA 

C4  = 

rg 

A7  5  = 

GHBHKr9K0^ 

BS  = 

C5  = 

0.0 

A76  = 

ghbhK0‘wh 

B6  = 

0.0 

C6  = 

^GA 

A7  7  - 

GHBHKr9K9ii;G/v  BH 

B7  = 

Ks 

Gh>ga 

C7  = 

0.0 

The  closed  loop  transfer  function  H(s)  relating  the 
missile  pitch  angle  9  to  the  commanded  pitch  angle  9c  is: 


H  (S)  = 


,y5  H.  N4s4  *  N3s3  +  N2i2  *  N^s  ♦  ,M0 _ 

s7  +  D.s6  +D.s5  +  D.s4  +  D-.s'5  +  D..S2  +  D,  s  +  D„ 
0  3  4  3  2  10 


where : 


N5  =  'M5+GhVGI!BH 


Z. 


N4  -  -  Vs  -M5BL*(Mc,T-M5T^GLBLBH*TTBLBH*TBLBH*GHBHT*GHBLBH 


n3  - 


N,  = 


tGLBLBHf*M<xGLBL*Mc.GHBH*GHBLBHTr*GHWL"BH 

n  ^  Z 

-M5WH2BL*1,H2'M,»T-M«-ir)-M6WL2BH*BLBH'MST-M 

*WL2  CMc,T-M5TT’  VGLBL^BHGLBlVBLGH  W  VGHBHTT 
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•  1 1 

i  1 


7  z 

"l  ■  -MS,*H2wL2»WH2BLtMa?'M5T^BHWL2tMc.T-M5:r)+WH2GLBLMa 


+W. “G„B„M 
L  H  H  a 


Z  Z 
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For  the  bending  mode  parameters  selected  and  the  plant 
parameters  for  time  =  40  seconds  the  transfer  function  H(s) 


t  * 


tI,„,_9(s)  _34 . 97  fs+  48 . 4)  (s-48.7)  (s+  0.018)  fs  +  0. 215  +  /-  j48 . 65) 

H  ( S) -Q-JSJ-  (-s  +  0.  o  17 (s' +  3 . "6T+  /  -  j  (r.~9 97 )( s~5\ '$'4 + j '3'8' . TS  j 

c  ( s -  2 . 91  +  / - j  62 . 25) 
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The  transfer  function  H(s)  was  factored  using  the  "roots" 
program  contained  in  Ref.  [7] .  The  location  of  the  poles  was 
confirmed  by  finding  the  eigenvalues  of  the  A  matrix  using 
the  "BASMAT"  program  also  contained  in  Ref.  [7].  The  poles 
and  zeros  of  H(s)  are  plotted  in  Figure  2-16. 

The  system  composed  of  the  plant  plus  bending  modes  is 
clearly  unstable  with  four  poles  of  the  closed  loop  transfer 
function  in  the  right  half  of  the  S-plane. 

Figure  2-17  shows  the  simulation  the  system  shown  in 
Figure  2-9  for  a  step  input  of  0.1  radians.  Line  1  of 
Figure  2-17  is  the  response  of  the  plant  without  bending 
modes  and  line  2  is  the  response  of  the  system  composed  of 
the  plant  with  bending  modes.  The  vertical  scaling  in  Figure 
2-17  was  restricted  to  keep  the  response  of  the  plant  in 
proper  perspective.  The  restriction  of  the  vertical  scale 
causes  saturation  of  the  plot  of  the  system  response  (line  2). 
Figure  2-18  shows  the  response  of  the  system  with  the  verti¬ 
cal  scale  extended  to  remove  the  effects  of  saturation  of 
the  plot.  For  clarity,  the  response  of  the  plant  without 
bending  modes  was  not  included  in  Figure  2-18. 

To  validate  the  s imulat ion , , the  time  response  of  the 
state  equations  with  a  0.1  radian  step  input  was  determined 
by  using  the  "RTRESP"  program  contained  in  Ref.  [7] .  The 
pitch  angle  of  the  missile,  9  was  found  to  be  described  by 
the  equation: 


Figure  2'16  Location  of  Poise  and  Zeros 
of  the  Plant  with  Two 
Bending  Modes 
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System;  Third  order  missile  model  with  two 
bending  modes  and  no  adaptive  notch 
filters 


Input:  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  two  bending 
modes 

Figure  2t17  Missile  Response  with  Two  Bending  Modes 
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System;  Third  order  missile  model  with  two 
bending  modes  and  no  adaptive  notch 
filters 

Input:  0.1  radian  step  command 


Figure  2-18  Missile  Response  with  Two  Bending 
Modes  (Expanded  Scale) 


_ _ 


0(t)  =  -0.09952  +  0.0457  e 

+  0.0582  e'3-63t  cos  (6 . 997t- 0 . 41 ) 

+  0.00246  eJ-a4t  cos (38 . 15t-2 . 45) 

+  0.00093  e2‘91t  cos (62 . 25t+3. 06) 

At  time  =  1.0  seconds,  the  above  equation  indicates  the 
pitch  angle  will  be  0(1)  =  -0.10162  radians.  The  numerical 
results  of  the  simulation  used  to  generate  Figure  2-17  indi¬ 
cate  the  pitch  angle  to  be  9(1)  =  -0.1017  radians.  The 
preceeding  results  indicate  the  simulation  program  is  valid. 

The  system  was  simulated  with  only  the  high  frequency 
bending  mode  and  with  only  the  low  frequency  bending  mode. 

The  purpose  of  these  simulations  was  to  show  that  each  bending 
mode  would  by  itself  destabilize  the  system  and  that  the  fre¬ 
quency  of  the  oscillation  of  the  resulting  unstable  system 
would  reflect  very  nearly  the  center  frequency  of  the 
bending  mode  that  destabilized  the  system.  Figure  2-19  illus¬ 
trates  the  response  of  the  system  for  a  0.1  radian  step 
input  when  only  the  high  frequency  bending  mode  was  present. 
Figure  2-20  shows  the  same  simulation  with  an  expanded  verti¬ 
cal  scale.  Figure  2-21  shows  the  response  of  the  system  for 
a  0.1  radian  step  input  with  only  the  low  frequency  bending 
mode  present.  Figure  2-22  shows  the  same  simulation  with 
an  expanded  vertical  scale.  Figure  2-23  shows  the  system 
response  with  both  bending  modes  for  a  ramp  input 
9  (t)  =  O.OSt. 
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System;  Third  order  missile  model  with  10  Hz 
bending  mode  and  no  adaptive  notch 
filter 


Input;  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  10  Hz  bending 
modes 

Figure  2-19  Missile  Response  with  10  Hz  Bending  Modes 


57 


Time,  Sec. 


System:  Third  order  missile  model  with  10  Hz 

bending  mode  and  no  adaptive  notch 
filter 

Input:  0.1  radian  step  command 


Figure  2-20  Missile  Response  with  10  Hz  Bending 
Mode  (expanded  scale) 
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System:  Third  order  missile  model  with  6  Hz 

bending  mode  and  no  adaptive  notch 
filters 

Input:  0,1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  6  Hz  bending 
modes 

Figure  2-21  Missile  Response  with  6  Hz  Bending  Mode 
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System;  Third  order  missile  model  with  6  H 
bending  mode  and  no  adaptive  notch 
filters 

Input:  0.1  radian  step  command 


Figure  2-22  Missile  Response  with  6  Hz  Bending 
Mode  (expanded  scale) 
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Figures  2-19  through  2-22  clearly  indicate  that  each 
bending  mode  is  in  itself  destabilizing.  Analysis  of  the 
numerical  data  associated  with  Figures  2-19  through  2-22 
confirms  that  the  frequency  of  the  oscillations  is  essen¬ 
tially  at  the  frequency  of  the  bending  mode  center  frequency. 


III.  ADAPTIVE  DIGITAL  FILTER  DESIGN 

A.  FACTORS  AFFECTING  THE  DESIGN 

The  emphasis  of  the  design  was  simplicity.  Significant 
improvement  in  the  results  could  have  been  achieved  by  use 
of  higher  order  filters  and  by  using  more  sophisticated  fil¬ 
ters,  but  the  improved  performance  would  have  been  accom¬ 
plished  at  the  expense  of  obscuring  the  very  simplicity  which 
makes  the  proposed  technique  so  appealing. 

The  sampling  frequency  used  throughout  this  thesis  is  100 
hertz.  This  selection  is  dictated  by  the  sampling  frequency 
used  by  the  autopilot  in  the  Trident  Missile.  Higher  sam¬ 
pling  frequencies  were  found  not  to  significantly  improve 
the  performance  of  the  algorithm,  however,  sampling  fre¬ 
quencies  lower  than  60  Hz  degraded  the  algorithm  to  where 
it  would  not  perform  reliably. 

The  bending  mode  frequencies  were  selected  such  that  one 
and  only  one  bending  mode  would  be  in  each  of  the  bands  of 
5-8  hertz  and  8-12  hertz  as  shown  in  Figure  2-6.  A  priori 
knowledge  of  the  band  or  range  of  frequencies  over  which  a 
given  bending  mode  may  vary  during  the  flight  of  the  missile 
is  essential  in  this  design.  If  specific  events  during  the 
flight  such  as  staging  change  the  bands  of  the  bending  modes 
then  the  algorithm  may  also  be  modified  to  adapt  to  the  new 
distribution  with  the  only  requirement  being  that  the  times 
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or  events  which  change  the  band  limitations  be  preprogrammed 
in  the  microprocessor  algorithm. 


B.  DESIGN  OF  THE  ALGORITHM  MODULES 

Figure  3-1  shows  the  modular  nature  of  the  algorithm  for 
the  determination,  tracking  and  notch  filtering  of  the  bending 
mode  signals.  The  sequence  of  modules  associated  with  the 
determination,  tracking  and  notch  filtering  for  each  of  the 
bending  modes  are  collectively  termed  a  channel.  The  channels 
are  identical  except  for  the  numerical  values  of  the  coeffi¬ 
cients;  therefore,  only  the  design  of  the  8-12  hertz  (high) 
channel  is  presented  in  detail.  The  results  of  the  design 
of  each  of  the  modules  comprising  the  low  channel  are  pre¬ 
sented  when  different  from  the  corresponding  module  in  the 
high  channel. 

1 .  The  Bandpass  Filter  Module 

The  overall  bandpass  filter  was  designed  as  two 
cascaded  second-order  bandpass  sections.  The  first  section 
is  a  first-order  Butterworth  lowpass  Filter  transformed  into 
a  bandpass  filter  and  then  transformed  into  the  z-domain 
using  the  bilinear  transformation.  The  second  section  is 
also  a  bandpass  filter  designed  to  adjust  the  overall  shape 
of  the  gain  curve  to  the  desired  form. 

A  first-order  Butterworth  lowpass  Filter  normalized 
to  a  cutoff  frequency  of  1  radian/second  has  a  transfer 
function : 
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Figure  3-1  Adaptive  Notch  Filter  Algorithm  Modules 


The  transformation  from  lowpass  to  bandpass  is: 


S  -  - Bs~^~ 

where : 

B  is  the  3  dB  bandwidth  defined  as  W  -  W, 

c  L 

W  is  the  center  frequency  defined  as  / 

W  is  the  upper  5  dB  frequency  in  rad/sec  and  WT  is  the 
lower  3  dB  frequency  in  rad/sec. 

The  resulting  analog  domain  transfer  function  for 
the  bandpass  filter  is: 


hbpCs)  - 


Bs 

- - - 

S  +  Bs  +  W  4 
c 


(4) 


Before  applying  the  bilinear  transformation  to  trans 
form  the  analog  bandpass  filter  to  the  discrete  domain  it 
is  necessary  to  prewarp  the  center  frequency  W  [8] .  The 
prewarping  is  done  by  using  the  following  transformation: 


WA  " 


“n 

2  TAN  (t) 


(5) 
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where : 


is  the  desired  discrete  or  digital  frequency  defined 

as : 

o  =  W  T 
“D  c 

T  is  the  sampling  frequency. 

1/2 

W  is  the  analog  center  frequency  defined  as  W  =  (W 

Selecting  8.0  hertz  as  the  lower  cutoff  frequency  and 
12.0  hertz  as  the  upper  cutoff  frequency.  The  values  to  be 
used  in  equations  (4)  and  (5)  become: 

=  2 tt -  8 . 0  =  50.265  rad/sec 

W  =  2-rr •  1 2 . 0  =  75.398  rad/sec 
u 

W  =  (W  W.)1/2  =  61.652  rad/sec 
c  u  L 

B  =  W  -W.  =  25.133  rad/sec 
u  L 

Applying  equation  (5)  with  the  selected  value  of  W 
yield  the  prewarped  center  frequency: 

=  63.582  rad/sec 

Substituting  the  above  value,  equation  (4)  becomes: 

H(s)  =  -2 - 2-5--1-5--  -  (6) 

s  +  25.133  s  +  4042.77 
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The  bilinear  transformation  for  the  s-domain  to  the 
z-domain  is: 


_  2 (Z  -  1) 
S  "  T(Z  +  1) 


Applying  the  transformation  to  equation  (6) ,  the  dis 
Crete  transfer  function  for  the  bandpass  filter  becomes: 


„  2 


H  (Z-)  =  0 . 103  (Z  -  1) 


Z  -1.4677  Z  +  0.79405 
After  slight  adjustments  to  obtain  the  desired  gain 
and  crossover  frequency  with  the  low  frequency  channel  band¬ 
pass  filter  (see  Figure  3-9),  the  transfer  function  H^(Z)  is 

0.1  (Z2  -  1) 


Vz)  ’  rr 


1.45625  Z  +  0.81 


The  pole  and  zero  configuration  of  H^CZ)  is  shown  in 
Figure  3-2  and  the  frequency  response  of  H^CZ)  is  shown  in 
Figure  3-3. 


As  expected,  the  frequency  response  of  H^(Z)  is  too 


flat  to  provide  the  needed  stopband  attenuation  for  the 
adjacent  channels.  The  frequency  response  of  the  overall 
bandpass  filter  is  greatly  improved  by  adding  a  second  sec¬ 
tion  of  filter  with  a  sharp  resonant  peak  at  the  center  fre¬ 
quency.  A  bandpass  filter  with  the  desired  resonant  peak 
was  designed  directly  in  the  z-domain  using  a  digital  filter 
frequency  response  program  (program  3-5)  to  adjust  the  poles 
and  zeros  to  obtain  the  desired  characteristics.  The  trans¬ 


fer  function  of  the  second  filter  section  is: 


„  _  0.456(1  -  0.809  Z  +  0.25) 

H  -  - 7 - - 


1.3753  Z  +  0.7225 
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Frequency,  Hz  (linear  scale) 


Figure  3- 


3  Frequency  Response  of  the  "irst  Section 
of  the  High  Frequency  Digital  Bandpass 
Filter 


The  pole  and  zero  configuration  of  H^CZ)  is  shown  in 
Figure  3-4.  The  frequency  response  of  H^CZ)  is  shown  in 
Figure  5-5. 

Figure  3-6  shows  the  flow  graph  of  the  algorithm  for 
the  implementation  of  the  two  cascaded  sections  and  Figure 
3-7  shows  the  resulting  frequency  response  of  the  cascaded 
sections . 

Figure  3-8  shows  the  frequency  response  of  the  low 
channel  bandpass  filter  designed  in  the  same  manner.  Figure 
3-9  illustrates  the  resulting  band  coverage  of  each  channel. 

The  possibility  of  a  bending  mode  occurring  at  the 
eight -hertz  cross-over  frequency  between  channels  causing 
both  channels  to  lock  on  that  frequency  and  leaving  one  bend¬ 
ing  mode  unattenuated  is  compensated  for  by  initializing  both 
notch  filters  at  the  cross-over  frequency,  thereby  attenuat¬ 
ing  a  bending  mode  at  eight-hertz  immediately  and  freeing 
either  channel  to  track  and  remove  the  unattenuated  bending 
mode . 

2 .  The  Zero  Crossing  Detector  Module 

The  zero  crossing  detector  is  based  on  a  linear 
estimate  of  the  crossing  time  if  the  amplitude  of  the  output 
of  the  bandpass  filter  changes  signs  between  samples.  Linear 
estimation  of  the  crossing  time  is  a  feature  incorporated 
into  the  digital  simulation  language  program  (DSL)  as  the 
cross  function  but  could  also  be  programmed  using  the 
algorithm  shown  in  Figure  3-10.  The  potential  for  using 


71 


N 


N 


r-H  O  O 

O  O  *“H 

H  II  >1 


lO 

CM 

vO 

O 


II 


oo 

o 

II 


iH  CM  rH  CM 
X  N  N  (N  N 

O  <  <  CQ  (tt 


G 

O 


c 

CD 

B 

o 

rH 

o, 

£ 


5-i 

(D 

•M 


PU 

rH 

CO 

■H 

•H 

00 

•H 

o 


O 

*• 

to 

0 


CTi 

to 

3 

O0 

LO  LO 

DO 

LO 

o 

cm 

*r-f 

\o 

00  LO 

to  CM 

u. 

•  CM 

• 

• 

o 

rH  • 

o 

i  o 

1  o 

(1 

II  II 

II  II 

rH  CM 

rH  CM 

*-4 

rH  rH 

rH  rH 

o 

A 

A 

B 

B 

Figure  3-10  Zero  Crossing  Time  Algorithm 


78 


ir 


a  sinusoidal  curve  estimation  rather  than  a  linear  estima¬ 
tion  was  considered  and  dismissed  because  any  increase  in 
the  accuracy  of  the  estimate  of  the  zero  crossing  time  was 
offset  by  the  jitter  in  the  actual  crossing  time  caused 
by  leakage  into  the  channel  of  the  adjacent  bending  mode  and 
plant  frequencies. 

The  major  deficiency  in  the  use  of  a  zero  crossing 
detector  as  a  frequency  measuring  device  is  the  nonlinear 
error  in  the  frequency  estimate  as  a  function  of  the  actual 
frequency.  For  example,  if  the  frequency  being  measured  is 
5.0  hertz  the  half-period  measured  should  be  0.1  seconds. 

If  noise  caused  a  0.01  second  measurement  error  and  the 
measured  half-period  was  0.09  seconds  then  the  frequency 
detector  would  estimate  the  input  frequency  to  be  5.55  hertz. 
If  the  input  frequency  was  10.0  hertz  the  same  0.01  second 
error  would  result  in  an  estimate  of  the  input  frequency  of 
12.5  hertz,  a  five  times  greater  error  for  the  same  measure¬ 
ment  error. 

The  effect  of  the  nonlinear  quality  of  the  estimated 
frequency  in  the  presence  of  noise  causes  a  much  greater 
variation  in  the  estimate  of  the  higher  frequency  bending 
modes  and  contributes  to  the  limit  cycle  problem  described 
below. 

3 .  The  Frequency  Averaging  Module 

The  primary  source  of  limit  cycle  behavior  in  the 
tracking  algorithm  appears  to  be  due  to  leakage  of  an 
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adjacent  bending  mode  frequency  into  the  zero  crossing  detec¬ 
tor  causing  a  periodic  error  in  the  frequency  determination. 
The  frequency  averaging  module  simply  determines  the  mean 
of  the  last  detected  frequency  from  the  frequency  detector 
and  the  current  notch  filter  center  frequency.  The  mean  of 
the  two  frequencies  becomes  the  new  notch  filter  center  fre¬ 
quency.  The  process  of  averaging  the  last  detected  frequency 
with  the  current  notch  filter  center  frequency  reduces  the 
amplitude  of  the  limit  cycle  but  does  not  eliminate  it.  To 
further  reduce  the  amplitude  of  the  limit  cycle,  several 
approaches  may  be  taken: 

(a)  Greatly  improve  the  stopband  attenuation  and  roll¬ 
off  characteristics  of  the  bandpass  filters.  This  approach 
is  discussed  in  the  conclusions  contained  in  Chapter  Five. 

(b)  Average  the  detected  frequency  over  a  sufficient 
number  of  zero  crossings  before  moving  the  notch.  This 
approach  was  found  to  significantly  reduce  the  amplitude  of 
the  limit  cycle  but  was  so  slow  in  moving  the  notch  to  the 
correct  frequency  that  the  bending  modes  would  increase  to 
such  an  amplitude  that  the  missile  would  be  destroyed 

4 .  Design  of  the  Notch  Filter  Module 

The  design  of  a  discrete  notch  filter  and  the  calcu¬ 
lation  of  the  coefficients  is  simple  and  efficient.  Figure 
3-11  illustrates  the  pole  and  zero  location  of  a  second- 
order  notch  filter  on  the  2-plane.  Since  the  zero  is 
located  on  the  unit  circle,  the  center  frequency  W  of  the 
notch  is  completely  specified  by  the  angle  n  . 
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FiL  TF|i  ROLE  ZERO  LOCATION 


Figure  3-11  General  Pole/Zero  Locations 
for  a  Second-Order  Discrete 
Notch  Filter 


The  angle  Qc  is  referred  to  as  the  discrete  fre¬ 
quency  and  is  defined  as  the  angle  from  the  positive  real 
axis  to  the  zero  such  that: 


ft  =  W  T 
c  c 


where : 


W  is  the  desired  notch  center  frequency  in  rad/sec 
T  is  the  sampling  period  in  seconds 

The  location  of  the  zero  on  the  unit  circle  at  a 


selected  discrete  frequency,  ft  ,  ensures  that  the  center 
of  the  notch  will  be  at  the  desired  frequency  completely 
independent  of  the  selected  location  of  the  poles.  The  loca¬ 
tion  of  the  poles  on  the  same  radial  (discrete  frequency) 
as  the  zero  maintains  the  approximate  symmetry  of  the  notch. 
The  radius  of  the  poles  determines  the  bandwidth  of  the 
notch.  The  closer  the  poles  are  to  the  unit  circle,  the 
narrower  the  notch  will  be  and  the  longer  the  settling  time. 

The  transfer  function  H^(Z)  for  the  notch  filter  is 
of  the  form: 


H  (Z)  =  -  A1  Z  +  1) 

N  Z2  -  B1  Z  +  B2 


where  the  coefficients  are  defined  by  the  following  formulas: 

A1  =  2  cos  (WcT) 

B!  =  2  r  cos  (W  T) 
p  c 

-> 

r:  =  r  " 

p 


r  is  the  radius  of  the  poles 
P 

W  is  the  desired  center  frequency  of  the  notch  in  rad/sec 
G  is  the  gain 

For  unity  DC  gain  H(Z)jZ=+l=l  yields  the  formula 
„  1  -  B1  +  B2 


The  relationship  between  the  pole  radius  and  the  3  dB 
bandwidth  of  the  notch  filter  is  a  function  of  center  fre¬ 
quency  but  can  be  approximated  by  making  certain  simplifying 
assumptions.  Figure  3-12  shows  the  geometry  on  the  unit 
circle  when  the  discrete  frequency,  W^T  is  at  the  3  dB 
cutoff  frequency.  The  gain  at  the  cutoff  frequency  is: 


r 

u 


0.707 


(7) 


In  this  application  it  is  reasonable  to  assume  that 
I Z 2 | ~ | P ^ | ,  the  poles  are  close  to  the  unit  circle,  and  the 
bandwidth  of  the  filter  is  small  (W  TCW^T) .  With  the  above 
assumptions,  equation  (7)  then  reduces  to: 


G 


0.707 


(8) 


2 

Using  the  Pythagorean  theorem  and  substituting  Zj  = 

2 

0.5  from  equation  (8)  the  following  relationship  can  be 
found : 


83 


(1  -  rp)  =  Z1 


(9) 


For  small  angles  the  cord  can  be  assumed  to  be  equal 
to  the  ARC  length,  therefore: 

Z1  =  WcT  -  WDT  «  O.S  BD  -  0.5  BWT  (10) 


where : 

Bp  is  the  discrete  bandwidth  in  rad 
Bw  is  the  analog  bandwidth  in  rad/sec 

Substituting  equation  (10)  into  equation  (9)  and 
solving  for  yields  the  approximate  relationship  between 
the  notch  filter  analog  frequency  bandwidth  and  the  pole 
radius  in  the  discrete  domain. 


BW  " 


2(1  - 


Ini 


(id 


The  numerical  data  associated  with  the  frequency 
response  curves  for  the  notch  filters  indicates  that  equation 
(11)  yields  bandwidths  that  are  larger  than  the  desired 
bandwidths  but  the  error  is  less  than  10  percent  in  the 
frequency  range  of  5  to  15  hertz. 

Figure  3-13  shows  the  gain  vs  frequency  response  of 
the  low  channel  notch  filter  with  a  center  frequency  of  6.0 
hertz  and  a  pole  radius  of  0.9.  This  configuration  has  a 
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3  dB  bandwidth  of  2.95  hertz.  Figure  3-14  shows  the  phase 
vs  frequency  response  of  the  same  notch  filter. 

Figure  3-15  shows  the  gain  vs  frequency  response  of 
the  high  channel  notch  filter  with  a  center  frequency  of 
10.0  hertz  and  a  pole  radius  of  0.8.  This  configuration  has 
a  3  dB  bandwidth  of  5.9  hertz.  Figure  3-16  shows  the  phase 
vs  frequency  response  of  the  same  notch. 

The  argument  used  to  create  Figure  2-6  established 
that  the  higher  of  the  two  bending  modes  would  have  a  wider 
bandwidth  than  the  lower  frequency  bending  mode,  therefore 
the  bandwidth  of  the  notch  filter  for  the  high  channel  was 
designed  to  be  wider  than  the  low  channel.  There  is  a  limit 
to  which  the  notch  bandwidth  can  be  increased  without  in¬ 
ducing  excessive  phase  loss  at  the  plant  operating  frequency 
of  two-hertz.  For  the  notch  bandwidths  used  in  this  algorithm 
the  high  channel,  with  the  notch  centered  at  10.0  hertz,  in¬ 
troduces  a  phase  loss  of  7.75  degrees  at  two-hertz  and  the 
low  channel,  with  the  notch  centered  at  6.0  hertz,  intro¬ 
duces  a  phase  loss  of  11.0  degrees  at  two-hertz.  The  com¬ 
bined  phase  loss  introduced  by  the  notch  filters  used  in 
this  algorithm  does  not  have  a  pronounced  effect  on  the 
stability  of  the  plant. 

5 .  Flow  Chart  of  the  Algorithm 

Figure  3-17  provides  a  flow  chart  of  the  entire 
algorithm.  The  "prime  time"  computations  are  the  processing 
of  the  error  signal  and  outputting  the  filtered  value  to  the 
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Figure  3-17  Algorithm  Flowchart 


91 


autopilot  with  minimum  processing  delay.  The  "slack  time" 
computations  are  the  computations  performed  after  the  error 
signal  has  been  processed  and  which  do  not  introduce  pro¬ 
cessing  delay  in  the  control  system.  With  a  100  hertz 
sampling  frequency,  sufficient  time  should  be  available  to 
complete  all  slack  time  calculations  every  cycle,  with 
currently  available  microprocessors. 

C.  TEST  OF  THE  ALGORITHM 

The  algorithm  was  first  tested  as  shown  in  Figure  3-18. 
An  input  signal  consisting  of  the  sum  of  equal  amplitude 
2.0,  6.0  and  10.0  hertz  sine  waves  initially  in  phase  was 
used  to  simulate  the  rate  gyro  output  signal  with  bending 
mode  contamination.  The  signal  was  first  passed  through 
the  tracking  notch  filter  modules  and  the  result  was  then 
used  as  an  input  to  the  two  channels  of  the  bending  mode 
frequency  determination,  tracking  and  filtering  algorithm. 
The  unusual  process  of  passing  the  input  signal  through  the 
notch  filters  first  was  done  to  allow  the  effects  of  feed¬ 
back  to  be  observed  as  the  notch  filters  remove  the  very 
signals  that  the  remainder  of  the  algorithm  is  trying  to 
track. 

Figure  3-19  shows  the  two-hertz  input  signal  (line  3) 
to  the  notch  filters  and  the  resulting  output  signal  (line 
2,  accentuated).  The  notch  filters  were  reasonably  able 
to  filter  out  the  6  and  10  hertz  components  and  pass  the 
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LOW  CHANNEL 
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Figure  3-18  Adaptive  Notch  Filter 

Algorithm  Test  Arrangement 


two-hertz.  The  numerical  data  associated  with  Figure  3-19 
indicates  an  average  phase  loss  of  18  degrees  was  experienced 
over  the  2  seconds  of  simulation. 

Figure  3-20  shows  the  tracking  history  of  each  of  the 
channels.  The  low  channel  locked  on  the  six-hertz  frequency 
(actual  lock  on  frequency  was  6.0546  Hz)  after  0.56  seconds 
and  remained  locked  on  for  the  remainder  of  the  simulation. 

The  high  channel  went  into  a  four  cycle  per  second  limit 
cycle  with  an  average  frequency  of  9.85  hertz.  The  cause 
of  the  limit  cycle  is  leadage  of  the  two-hertz  signal  when 
the  amplitude  of  the  two-hertz  reaches  near  maximum  ampli¬ 
tudes.  The  limit  cycle  has  a  minimum  effect  on  the  output 
of  the  notch  filters  since  the  high  channel  notch  bandwidth 
is  nearly  six-hertz  and  the  magnitude  of  the  limit  cycle 
is  only  1.65  hertz. 

Figure  3-21  shows  the  input  (line  1)  and  output  (line  2, 
accentuated)  of  the  low  frequency  channel  bandpass  filter. 
Figure  3-22  shows  the  input  (line  1)  and  the  output  (line  2, 
accentuated)  of  the  high  channel  bandpass  filter.  The  effects 
of  the  leadage  of  the  two-hertz  signal  can  be  seen  in  the 
output  of  both  filters. 
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Figure  3-20  Frequency  Detector  Output  for  6  and 
10  Hz  Input 
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Line  1  is  the  output  from  the  adaptive  notch 
filter 


Line  2  is  the  output  from  the  low  channel 
bandpass  filter 


Figure  3-21  Low  Channel  Bandpass  Filter  Output 


Time,  Sec. 

Line  1  is  the  output  from  the  adaptive  notch 
filters 

Line  2  is  the  output  from  the  high  channel 
bandpass  filter 


Figure  3-22  Output  from  the  High  Channel  Bandpass 
Filter 


IV.  SIMULATION  STUDIES 


A.  DIGITAL  SIMULATION  LANGUAGE  (DSL/360) 

The  programming  languages  used  throughout  this  thesis  are 
the  IBM  Digital  Simulation  Language  (DSL/360)  [9]  and  FORTRAN 
IV.  DSL/360  has  not  been  released  for  general  use  by  the  IBM 
Corporation  but  was  used  in  this  thesis  for  the  simulation 
studies  due  to  its  enhanced  capabilities.  The  exclusive 
features  of  DSL/360  were  used  only  in  modeling  of  the  con¬ 
tinuous  systems.  The  modules  of  the  adaptive  digital  filter 
are  programmed  exclusively  in  FORTRAN  IV  and  may  be  used  with 
any  continuous  systems  simulation  language  that  will  accept 
FORTRAN  subroutines. 

B.  BLOCK  DIAGRAM  OF  THE  SIMULATION  MODEL 

Figure  4-1  illustrates  the  block  diagram  of  the  system 
used  to  develop  Program  4-2.  The  block  diagram  forms  the 
core  program  for  all  the  simulation  studies  in  this  chapter 
with  the  only  changes  being  the  addition  of  a  ramp  function 
for  the  input  command  in  Program  4-3  and  ramp  functions  to 
create  dynamic  bending  mode  frequencies  in  Programs  4-10  and 
4-12. 

C.  SIMULATION  CONDITIONS 

The  parameters  of  the  plant  selected  for  all  the  following 
simulation  studies  are  from  Table  2-1  at  time  =  40  seconds. 
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The  selected  plant  parameters  represent  the  values  when  the 
missile  passes  through  the  point  of  maximum  dynamic  pres¬ 
sure  (max  Q)  and  is  the  least  stable  even  without  the  pre¬ 
sence  of  bending  modes.  When  simulation  studies  were 
conducted  with  plant  parameters  yielding  a  more  stable  plant 
the  results  indicated  the  proposed  adaptive  digital  filter¬ 
ing  technique  functioned  essentially  the  same  and  therefore 
only  the  worst  case  conditions  are  presented. 

Two  pitch  commands  from  the  guidance  system  were  selected 
for  simulation:  a  0.1  radian  step  input  and  a  0.05  radian 
per  second  ramp.  The  step  command  simulates  a  command  from 
the  autopilot  which  would  normally  occur  at  water  exit, 
staging,  or  inflight  restart  of  the  guidance  system  after  a 
failure.  The  simulation  studies  with  a  step  input  are  in¬ 
tended  to  simulate  the  condition  of  a  guidance  system 
restart  at  the  point  in  the  flight  when  the  missile  has  the 
least  stability,  a  highly  unlikely  but  undoubtably  severe 
condition.  The  ramp  input  simulates  a  more  normal  condition 
where  the  autopilot  responds  to  wind  shear  at  the  point  of 
least  stability.  The  ramp  input  response  serves  to  demon¬ 
strate  that  although  the  bending  modes  are  not  excited  as 
violently  as  they  are  with  a  step  input,  the  adaptive  digital 
filters  successfully  detect  their  presence  and  provide  the 
proper  filtering.  The  ramp  input  again  includes  the  simu¬ 
lation  of  an  inflight  restart  of  the  guidance  system  since 
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the  simulation  is  started  with  no  history  available  to  the 
adaptive  digital  filters  or  frequency  detector  algorithms. 

The  bending  mode  frequencies  were  the  primary  variable 
in  the  simulation  studies.  The  first  two  studies  have  fixed 
bending  mode  frequencies  centered  at  6.0  and  10.0  hertz. 

The  third  study  has  the  bending  mode  frequencies  again  fixed 
but  at  7.0  and  9.0  hertz  to  ensure  that  the  algorithm  can 
separate  two  close  bending  mode  frequencies.  The  fourth 
study  has  the  bending  modes  fixed  at  5.5  and  14.0  hertz  to 
test  the  ability  of  the  algorithm  to  properly  filter  two 
widely  separated  bending  mode  frequencies.  The  final  two 
studies  have  dynamic  bending  modes  with  the  upper  bending 
mode  changing  at  the  rate  of  1.5  hertz  per  second  and  the 
lower  bending  mode  changing  at  1.0  hertz  per  second.  The 
rates  of  change  far  exceed  the  rates  of  change  anticipated 
but  were  included  as  a  demonstration  of  the  capability  of 
the  algorithm  to  track  very  high  rates  of  change  in  the 
bending  mode  frequencies.  The  final  two  studies  illustrate 
the  response  of  the  system  with  the  bending  mode  frequencies 
converging  and  diverging. 

Table  4-1  provides  a  summary  of  the  conditions  used  in 
the  six  simulation  studies. 

D.  SIMULATION  RESULTS 

Figures  4-2  and  4-3  illustrate  the  system  response  for  a 
0.1  radian  step  pitch  command  with  the  bending  mode  fre¬ 
quencies  fixed  at  6.0  and  10.0  Hz.  Line  1  in  Figure  4-2  is 


102 


TABLE  4-1 


FIGURES 

INPUT 

BENDING  MODES 

4-2 

0.1  rad 

6.0  Hz  and  10.0  Hz 

4-3 

STEP 

FIXED 

4-4 

0.05*T  rad/sec 

6.0  Hz  and  10.0  Hz 

4-5 

RAMP 

FIXED 

4-6 

0.1  rad 

7.0  Hz  and  9.0  Hz 

4-7 

STEP 

FIXED 

4-8 

0 . 1  rad 

5.5  Hz  and  14.0  Hz 

4-9 

STEP 

FIXED 

4-10 

0.1  rad 

6.0  Hz  +  1.0  Hz/sec 

4-11 

STEP 

DYNAMIC 

4-12 

0.1  rad 

8.0  Hz  -  1.0  Hz/sec 

4-13 

STEP 

9.0  Hz  +  1.5  Hz/sec 

DYNAMIC 
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Pitch 

Angle 


'0. 00  O.iO  0.  ,0  0  SO  0.80  1-00  1.2Q  l  UQ  I  SO  1 .80  2.00 

Time,  Sec. 

System:  Third  order  missile  model  with  fixed 

bending  modes  at  6  and  10  Hz 


Input:  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 

Figure  4-2  Missile  Step  Response  with  Fixed  6  and 
10  Hz  Bending  Modes 
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Figure  4-3  Frequency  Detector  Output  for  Fixed 
6  and  10  Hz  Bending  Modes 


the  nominal  response  of  the  plant  without  bending  inodes  to 
the  same  input.  It  is  provided  as  a  standard  by  which  the 
effectiveness  of  the  adaptive  digital  notch  filtering  may 
be  evaluated.  Line  2  in  Figure  4-2  is  the  response  of  the 
system  with  the  bending  modes  and  the  adaptive  digital  notch 
filters . 

Figure  4-3  provides  a  history  of  the  center  frequency  of 
the  high  and  low  channel  notch  filters.  Line  1  is  the  his¬ 
tory  of  the  low  channel  [5. 0-8.0  He)  notch  center  frequency 
and  Line  2  is  the  history  of  the  high  channel  (8.0-12.0  Hz) 
notch  center  frequency.  Lines  3  and  4  indicate  the  actual 
bending  mode  center  frequencies. 

Figure  4-2  must  be  compared  with  Figure  2-17  to  gain  the 
proper  perspective  on  the  degree  to  which  the  system  has  been 
stabilized  by  the  adaptive  notch  filters.  Figure  2-17  is 
precisely  the  same  system  and  the  same  input  except  the 
adaptive  digital  notch  filters  were  not  in  the  system  simu¬ 
lated  in  Figure  2-17.  The  amount  of  leakage  of  the  6.0 
hertz  bending  mode  frequency  through  the  notch  filer  is  due 
to  the  relatively  low  precision  coefficients  used  in  the 
notch  filters.  Simulation  with  double  precision  coefficients 
reduced  the  leakage  considerably  but  the  implementation  of 
the  notch  filter  with  this  precision  was  not  considered 
practical  for  implementation  with  a  microprocessor  and 
therefore  not  pursued. 
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Figure  4-3  illustrates  that  the  algorithm  can  rapidly  and 


accurately  track  the  6.0  hertz  bending  mode  frequency  but 
exhibits  limit  cycle  behavior  when  tracking  the  10.0  hertz 
bending  mode.  The  limit  cycle  behavior  can  be  attributed  to 
the  low  frequencies  associated  with  the  plant  dynamics  causing 
errors  in  the  frequency  detector  as  discussed  in  Chapter 
Three.  The  limit  cycle  does  not  appear  to  diminish  the  high 
channel's  effectiveness  in  eliminating  the  10.0  Hz  bending 
mode  in  that  there  is  no  evidence  of  the  10.0  hertz  bending 
mode  in  the  output  and  it  is  attenuated  to  the  point  that  the 
frequency  detector  can  no  longer  track  it  after  1.5  seconds. 

Figures  4-4  and  4-5  show  the  response  of  the  system  for 
a  0.05*T  ramp  pitch  command  with  bending  modes  centered  at 
6.0  and  10.0  hertz.  Comparing  Figure  4-4  with  Figure  2-23 
it  is  evident  that  when  the  adaptive  digital  notch  filters 
are  included  in  the  system  the  response  is  essentially  the 
same  as  if  no  bending  modes  were  present.  Figure  4-5  shows 
the  ability  of  the  high  channel  to  much  more  accurately  lock 
onto  the  10.0  hertz  bending  mode  when  the  high  plant  rates 
are  not  present.  The  very  satisfactory  performance  for  a 
ramp  input  indicated  that  this  input  did  not  stress  the 
ability  of  the  adaptive  digital  notch  filter  algorithm  and 
it  was  therefore  not  used  in  the  remainder  of  the  simulation 
studies . 

Figures  4-6  and  4-7  show  the  response  of  the  system  for 
a  0.1  radian  step  input  when  the  bending  mode  frequencies 
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Time,  Sec. 

System:  Third  order  missile  model  with  fixed 

bending  modes  at  6  and  10  Hz 

Input:  0.05  radian/sec  ramp  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 

Figure  4^4  Missile  Ramp  Response  with  Fixed  6  and  10  Hz 
Bending  Modes 
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Figure  4-5  Frequency  Detector  Output  for  Fixed 
6  and  10  Hz  Bending  Modes 
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System:  Third  order  missile  model  with  fixed 

bending  modes  at  7  and  9  Hz 

Input:  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 


Figure  4-6  Missile  Step  Response  for  Fixed  7  and 
9  Hz  Bending  Modes 
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Figure  4-7  Frequency  Detector  Output  for  Fixed 
7  and  9  Hz  Bending  Modes 
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were  fixed  and  only  separated  by  2.0  hertz.  The  system 
response  is  very  close  to  the  nominal  plant  with  some  evi¬ 
dence  of  both  bending  modes  present  in  the  output.  The 
cross  channel  leakage  can  be  seen  to  cause  oscillations  in 
both  of  the  frequency  detectors. 

Figures  4-8  and  4-9  show  the  response  of  the  system  to 
a  0.1  radian  step  input  with  the  bending  modes  separated  by 
7.5  hertz.  The  lower  frequency  bending  mode  is  at  5.5  hertz 
or  0.5  hertz  above  the  lower  limit  to  which  the  algorithm 
is  allowed  to  track  and  the  upper  frequency  bending  mode  is 
14.0  hertz.  The  upper  bending  mode  frequency  is  outside  of 
the  passband  for  the  high  channel  frequency  detector  but 
within  the  15.0  hertz  that  the  high  channel  is  allowed  to 
track  a  bending  mode.  The  performance  of  the  system  is 
marginal  at  best.  This  simulation  clearly  illustrates  the 
inability  of  the  algorithm  to  reliably  acquire  or  track  a 
bending  mode  higher  than  12.0  hertz.  The  problem  is  not 
accentuated  by  the  bandpass  filter  in  the  high  channel 
attenuating  frequencies  above  12.0  hertz.  Increasing  the 
bandwidth  of  the  high  channel  to  15.0  hertz  did  not  enhance 
the  ability  of  the  algorithm  to  reliably  identify  a  fre¬ 
quency  above  12.0  hertz. 

Figures  4-10  and  4-11  show  the  system  response  for  a  0.1 
radian  step  input  with  dynamic,  converging  bending  mode 
frequencies.  The  high  frequency  bending  mode  starts  at  12.0 
hertz  and  decreases  at  1.5  hertz  per  second  to  9.0  hertz. 
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System:  Third  order  missile  model  with  fixed 

bending  modes  at  5.5  and  14.0  Hz 

Input:  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 

4-8  Missile  Response  with  Fixed  5.5  and  14.0  Hz 
Bending  Modes 
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ACTUAL  BENDING  MODE  FREQUENCY 


Frequency 
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Figure  4-9  Frequency  Detector  Output  for  Fixed 
5.5  and  14.0  Hz  Bending  Modes 
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System;  Third  order  missile  mode  with  converging 
dynamic  bending  modes 

Input;  0.1  radian  step  command 

Line  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 

Figure  4-10  Missile  Step  Response  with  Converging 
Dynamic  Bending  Modes 
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Figure  4-11  Frequency  Detector  Output  for  Converging 
Dynamic  Bending  Modes 
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The  low  frequency  bending  mode  starts  at  6.0  hertz  and  in¬ 
creases  at  1.0  hertz  per  second  to  8.0  hertz.  The  ability 
of  the  adaptive  digital  notch  filters  to  stabilize  the  system 
with  such  high  rates  of  change  was  not  expected.  When 
higher  rates  of  change  were  simulated  for  the  converging 
case,  the  error  in  the  high  frequency  estimate  became  so  poor 
that  the  system  became  unstable  at  the  frequency  of  the  high 
bending  mode.  A  rate  of  1.5  hertz  per  second  is  felt  to  be 
the  limit  on  the  rate  of  change  of  either  bending  mode  fre¬ 
quency. 

Figures  4-12  and  4-13  show  the  system  response  to  a  0.1 
radian  step  input  with  dynamic  diverging  bending  mode  fre¬ 
quencies.  The  upper  bending  mode  frequency  starts  at  9.0 
hertz  and  increases  at  the  rate  of  1.5  hertz  per  second  to 
12.0  hertz.  The  lower  frequency  bending  mode  starts  at  8.0 
hertz  and  decreases  at  the  rate  of  1.0  hertz  per  second  to 
6,0  hertz.  The  very  good  performance  of  the  system  is  pro¬ 
bably  due  to  the  notch  filters  being  initialized  at  8.0 
hertz.  The  inability  of  the  high  channel  to  track  a  bending 
mode  frequency  greater  than  12.0  hertz  is  evident  as  the 
high  bending  mode  frequency  increases. 
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Pitch 
Angle 
0 ,  Rad. 


r  iv  O.MO  0.60  U  .80  1  JC  WO  1 *0  TTsO  T.8<J  7. 00 

Time,  Sec. 

System:  Third  order  missile  model  with  diverging 

dynamic  bending  modes 

Input:  0,1  radian  step  command 

Lin-e  1  is  the  nominal  missile  response  without 
bending  modes 

Line  2  is  the  missile  response  with  bending  modes 


Figure  4-12  Missile  Step  Response  with  Diverging 
Dynamic  Bending  Modes 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS  * 

The  feasibility  of  implementing  an  adaptive  digital  notch 
filter  to  stabilize  the  Trident  Missile  Control  System  when 
two  dynamic  bending  modes  are  present  has  been  demonstrated. 
The  conditions  simulated  in  Chapter  Four  represent  the  worst 
case  conditions  and  the  performance  of  the  algorithm  under 
these  conditions  was  satisfactory  considering  the  very  de¬ 
sirable  response  to  less  violent  and  more  likely  commands 
as  shown  in  Figure  4-4. 

The  algorithms  used  in  the  simulation  studies  contained 
in  Chapter  Four  were  developed  with  the  ability  to  implement 
them  in  a  microprocessor  program  as  a  primary  concern.  Sim¬ 
plicity  of  the  algorithm  was  maintained,  perhaps  at  the 
expense  of  performance,  to  ensure  that  the  algorithm  could 
be  implemented  on  a  microprocessor  for  subsequent  real-time 
studies.  The  time  varying,  nonlinear  characteristics  of  the 
system  when  the  adaptive  digital  notch  filters  are  tracking 
dynamic  bending  modes  with  varying  input  commands  deny  an 
analytical  approach  and  therefore  require  extensive  simula¬ 
tion  to  optimize  the  algorithm. 
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B.  RECOMMENDATIONS  FOR  FURTHER  STUDIES 

The  simplicity  of  the  digital  bandpass  filters  allows  an 
excessive  amount  of  out  of  band  interference.  The  perfor¬ 
mance  of  the  bandpass  filters  could  be  significantly  improved 
without  an  appreciable  increase  in  the  order  of  the  filters 
by  using  more  sophisticated  types  of  filters  such  as  a 
Chebyshev  or  an  elliptical  filter.  The  temptation  to  enhance 
the  performance  of  the  bandpass  filter  in  nonreal-time  simu¬ 
lations  without  considering  the  impl  mentation  problem  on  a 
microprocessor  must  be  resisted. 

The  direct  realization  of  the  bandpass  and  notch  digital 
filters  did  not  minimize  memory,  noise  due  to  finite  pre¬ 
cision  coefficients  or  the  number  of  multiplications  required 
for  each  output.  The  structure  of  the  filters  needs  to  be 
examined  with  a  specific  microprocessor’s  capabilities  con¬ 
sidered.  The  2IL0G  Z-8000  microprocessor  capabilities  were 
considered  during  the  development  of  the  algorithm  as  a  cur¬ 
rent  state-of-the-art  microprocessor. 

The  algorithm  must  be  evaluated  on  the  six  degree  of 
freedom  model  of  the  Trident  Missile  available  at  Lockheed 
Missile  and  Space  Company,  Sunnyvale,  California.  The 
results  of  the  simulation  studies  contained  in  Chapter  Four 
are  only  as  valid  as  the  assumptions  made  in  modeling  the 
bending  modes  and  in  simplifying  the  plant.  It  must  be 
emphasized  that  the  algorithm  contained  in  this  thesis  will 
not  perform  satisfactorily  with  the  bending  mode  frequencies 
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present  in  the  current  Trident  Missile.  The  algorithm  in 
this  thesis  was  designed  to  identify,  track  and  attenuate 
bending  modes  that  are  a  full  octave  below  those  of  the 
current  Trident  Missile. 

A  phaselock  loop  scheme  should  be  reconsidered  as  a 
frequency  detector.  The  phaselock  loop  has  precisely  the 
desired  characteristic  of  seeking  out  and  locking  onto  a 
frequency  within  its  predetermined  acquisition  range  and 
seems  ideally  suited  for  the  task  of  i  mtifying  and  track¬ 
ing  the  bending  mode  frequencies.  It  is  recommended  that 
future  studies  of  the  application  of  phaselock  loops  in 
this  application  use  actual  phaselock  loops  rather  than  com¬ 
puter  simulation  due  to  the  inability  to  properly  simulate 
the  non  linear  behavior  of  the  phaselock  loop  with  algorithms 
known  to  the  author.  The  feasibility  of  prefiltering  the 
input  to  the  phaselock  loop  with  a  bandpass  filter  to  atten¬ 
uate  out  of  track  range  signals  and  thereby  reduce  the  lock- 
on  time  should  be  considered. 
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INITIAL  DIGITAL  BENDING  MODE  FILTER  PARAMETERS 
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MODEL  OF  THE  THIRD  ORDER  PLANT 
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INITIAL  DIGITAL  BENDING  MODE  FILTER  PARAMETERS 
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HZ  DIGITAL  BANDPASS  FILTER  AND  FREQUENCY  DETECTOR 
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COMPUTE  THE  CENTER  FREQUENCY  FOR  THE  HIGH  FREQUENCY  NOTCH 
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MODEL  OF  THE  THIRD  ORDER  PLANT 
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ASSOCIATED  FIGURES 
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HZ  DIGITAL  BANDPASS  FILTER  AND  FREQUENCY  QETECTOR 
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COMPUTE  THE  CENTER  FREQUENCY  FOR  THE  HIGH  FREQUENCY  NOTCH 
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MODEL  QF  THE  PLANT  WITHOUT  BENDING  MODES 
TO  ESTABLISH  NOMINAL  PERFORMANCE  . 
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ASSOCIATED  FIGURES 
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